using System;
using System.Collections.Generic;
using System.Text;
using NUnit.Framework;
using dnAnalytics.Math;
using dnAnalytics.LinearAlgebra;
using dnAnalytics.LinearAlgebra.Decomposition;
using dnAnalytics.Exceptions;

namespace dnAnalytics.UnitTests.LinearAlgebra.Decomposition
{
    [TestFixture]
    [Category("Managed")]
    public class LevinsonTest
    {
        #region Fields

        // unit testing - order 1

        Vector TR1;			// Toeplitz matrix
        Vector LC1;			// Toeplitz matrix
        Matrix A1;			// Lower triangle matrix
        Vector D1;			// diagonal vector
        Matrix B1;			// upper triangle matrix
        Matrix I1;			// inverse matrix
        double Det1;		// exact determinant
        Vector X1;			// RHS vector
        Vector Y1;			// LHS vector
        double Tolerance1;	// allowable tolerance

        // unit testing - order 2

        Vector TR2;			// Toeplitz matrix
        Vector LC2;			// Toeplitz matrix
        Matrix A2;			// Lower triangle matrix
        Vector D2;			// diagonal vector
        Matrix B2;			// upper triangle matrix
        Matrix I2;			// inverse matrix
        double Det2;		// exact determinant
        Vector X2;			// RHS vector
        Vector Y2;			// LHS vector
        double Tolerance2;	// allowable tolerance

        // unit testing - order 3

        Vector TR3;			// Toeplitz matrix
        Vector LC3;			// Toeplitz matrix
        Matrix A3;			// Lower triangle matrix
        Vector D3;			// diagonal vector
        Matrix B3;			// upper triangle matrix
        Matrix I3;			// inverse matrix
        double Det3;		// exact determinant
        Vector X3;			// RHS vector
        Vector Y3;			// LHS vector
        double Tolerance3;	// allowable tolerance

        // unit testing - order 4

        Vector TR4;			// Toeplitz matrix
        Vector LC4;			// Toeplitz matrix
        Matrix A4;			// Lower triangle matrix
        Vector D4;			// diagonal vector
        Matrix B4;			// upper triangle matrix
        Matrix I4;			// inverse matrix
        double Det4;		// exact determinant
        Vector X4;			// RHS vector
        Vector Y4;			// LHS vector
        double Tolerance4;	// allowable tolerance

        // unit testing - order 5

        Vector TR5;			// Toeplitz matrix
        Vector LC5;			// Toeplitz matrix
        Matrix A5;			// Lower triangle matrix
        Vector D5;			// diagonal vector
        Matrix B5;			// upper triangle matrix
        Matrix I5;			// inverse matrix
        double Det5;		// exact determinant
        Vector X5;			// RHS vector
        Vector Y5;			// LHS vector
        double Tolerance5;	// allowable tolerance

        // unit testing - order 10

        Vector TR10;		// Toeplitz matrix
        Vector LC10;		// Toeplitz matrix
        Matrix A10;			// Lower triangle matrix
        Vector D10;			// diagonal vector
        Matrix B10;			// upper triangle matrix
        Matrix I10;			// inverse matrix
        double Det10;		// exact determinant
        Vector X10;			// RHS vector
        Vector Y10;			// LHS vector
        double Tolerance10;	// allowable tolerance

        #endregion Fields

        #region Test Fixture Setup

        [TestFixtureSetUp]
        public void SetupTestCases()
        {
            // unit testing values - order 1

            TR1 = VectorBuilder.CreateVector(1, VectorType.Dense);
            TR1[0] = +3.0000000000000000E+000;

            LC1 = VectorBuilder.CreateVector(1, VectorType.Dense);
            LC1[0] = +3.0000000000000000E+000;

            A1 = MatrixBuilder.CreateMatrix(1, MatrixType.Dense);
            A1[0, 0] = +1.0000000000000000E+000;

            D1 = VectorBuilder.CreateVector(1, VectorType.Dense);
            D1[0] = +3.3333333333333331E-001;

            B1 = MatrixBuilder.CreateMatrix(1, MatrixType.Dense);
            B1[0, 0] = +1.0000000000000000E+000;

            Det1 = +3.0000000000000000E+000;

            I1 = MatrixBuilder.CreateMatrix(1, MatrixType.Dense);
            I1[0, 0] = +3.3333333333333331E-001;

            X1 = VectorBuilder.CreateVector(1, VectorType.Dense);
            X1[0] = +1.0000000000000000E+000;

            Y1 = VectorBuilder.CreateVector(1, VectorType.Dense);
            Y1[0] = +3.0000000000000000E+000;


            // unit testing values - order 2

            TR2 = VectorBuilder.CreateVector(2, VectorType.Dense);
            TR2[0] = +3.0000000000000000E+000;
            TR2[1] = +2.0000000000000000E+000;

            LC2 = VectorBuilder.CreateVector(2, VectorType.Dense);
            LC2[0] = +3.0000000000000000E+000;
            LC2[1] = +1.0000000000000000E+000;

            A2 = MatrixBuilder.CreateMatrix(2, MatrixType.Dense);
            A2[0, 0] = +1.0000000000000000E+000;
            A2[1, 0] = -3.3333333333333331E-001;
            A2[1, 1] = +1.0000000000000000E+000;

            D2 = VectorBuilder.CreateVector(2, VectorType.Dense);
            D2[0] = +3.3333333333333331E-001;
            D2[1] = +4.2857142857142855E-001;

            B2 = MatrixBuilder.CreateMatrix(2, MatrixType.Dense);
            B2[0, 0] = +1.0000000000000000E+000;
            B2[0, 1] = -6.6666666666666663E-001;
            B2[1, 1] = +1.0000000000000000E+000;

            Det2 = +7.0000000000000000E+000;

            I2 = MatrixBuilder.CreateMatrix(2, MatrixType.Dense);
            I2[0, 0] = +4.2857142857142855E-001;
            I2[0, 1] = -2.8571428571428570E-001;
            I2[1, 0] = -1.4285714285714285E-001;
            I2[1, 1] = +4.2857142857142855E-001;

            X2 = VectorBuilder.CreateVector(2, VectorType.Dense);
            X2[0] = +1.0000000000000000E+000;
            X2[1] = +2.0000000000000000E+000;

            Y2 = VectorBuilder.CreateVector(2, VectorType.Dense);
            Y2[0] = +7.0000000000000000E+000;
            Y2[1] = +7.0000000000000000E+000;


            // unit testing values - order 3

            TR3 = VectorBuilder.CreateVector(3, VectorType.Dense);
            TR3[0] = +3.0000000000000000E+000;
            TR3[1] = +2.0000000000000000E+000;
            TR3[2] = +1.0000000000000000E+000;

            LC3 = VectorBuilder.CreateVector(3, VectorType.Dense);
            LC3[0] = +3.0000000000000000E+000;
            LC3[1] = +1.0000000000000000E+000;
            LC3[2] = +0.0000000000000000E+000;

            A3 = MatrixBuilder.CreateMatrix(3, MatrixType.Dense);
            A3[0, 0] = +1.0000000000000000E+000;
            A3[1, 0] = -3.3333333333333331E-001;
            A3[1, 1] = +1.0000000000000000E+000;
            A3[2, 0] = +1.4285714285714285E-001;
            A3[2, 1] = -4.2857142857142855E-001;
            A3[2, 2] = +1.0000000000000000E+000;

            D3 = VectorBuilder.CreateVector(3, VectorType.Dense);
            D3[0] = +3.3333333333333331E-001;
            D3[1] = +4.2857142857142855E-001;
            D3[2] = +4.3750000000000000E-001;

            B3 = MatrixBuilder.CreateMatrix(3, MatrixType.Dense);
            B3[0, 0] = +1.0000000000000000E+000;
            B3[0, 1] = -6.6666666666666663E-001;
            B3[1, 1] = +1.0000000000000000E+000;
            B3[0, 2] = +1.4285714285714285E-001;
            B3[1, 2] = -7.1428571428571430E-001;
            B3[2, 2] = +1.0000000000000000E+000;

            Det3 = +1.6000000000000000E+001;

            I3 = MatrixBuilder.CreateMatrix(3, MatrixType.Dense);
            I3[0, 0] = +4.3750000000000000E-001;
            I3[0, 1] = -3.1250000000000000E-001;
            I3[0, 2] = +6.2500000000000000E-002;
            I3[1, 0] = -1.8750000000000000E-001;
            I3[1, 1] = +5.6250000000000000E-001;
            I3[1, 2] = -3.1250000000000000E-001;
            I3[2, 0] = +6.2500000000000000E-002;
            I3[2, 1] = -1.8750000000000000E-001;
            I3[2, 2] = +4.3750000000000000E-001;

            X3 = VectorBuilder.CreateVector(3, VectorType.Dense);
            X3[0] = +1.0000000000000000E+000;
            X3[1] = +2.0000000000000000E+000;
            X3[2] = +3.0000000000000000E+000;

            Y3 = VectorBuilder.CreateVector(3, VectorType.Dense);
            Y3[0] = +1.0000000000000000E+001;
            Y3[1] = +1.3000000000000000E+001;
            Y3[2] = +1.1000000000000000E+001;


            // unit testing values - order 4

            TR4 = VectorBuilder.CreateVector(4, VectorType.Dense);
            TR4[0] = +4.0000000000000000E+000;
            TR4[1] = +3.0000000000000000E+000;
            TR4[2] = +2.0000000000000000E+000;
            TR4[3] = +1.0000000000000000E+000;

            LC4 = VectorBuilder.CreateVector(4, VectorType.Dense);
            LC4[0] = +4.0000000000000000E+000;
            LC4[1] = +1.0000000000000000E+000;
            LC4[2] = +2.0000000000000000E+000;
            LC4[3] = +3.0000000000000000E+000;

            A4 = MatrixBuilder.CreateMatrix(4, MatrixType.Dense);
            A4[0, 0] = +1.0000000000000000E+000;
            A4[1, 0] = -2.5000000000000000E-001;
            A4[1, 1] = +1.0000000000000000E+000;
            A4[2, 0] = -5.3846153846153844E-001;
            A4[2, 1] = +1.5384615384615385E-001;
            A4[2, 2] = +1.0000000000000000E+000;
            A4[3, 0] = -8.1818181818181823E-001;
            A4[3, 1] = +9.0909090909090912E-002;
            A4[3, 2] = +9.0909090909090912E-002;
            A4[3, 3] = +1.0000000000000000E+000;

            D4 = VectorBuilder.CreateVector(4, VectorType.Dense);
            D4[0] = +2.5000000000000000E-001;
            D4[1] = +3.0769230769230771E-001;
            D4[2] = +2.9545454545454547E-001;
            D4[3] = +2.7500000000000002E-001;

            B4 = MatrixBuilder.CreateMatrix(4, MatrixType.Dense);
            B4[0, 0] = +1.0000000000000000E+000;
            B4[0, 1] = -7.5000000000000000E-001;
            B4[1, 1] = +1.0000000000000000E+000;
            B4[0, 2] = +7.6923076923076927E-002;
            B4[1, 2] = -7.6923076923076927E-001;
            B4[2, 2] = +1.0000000000000000E+000;
            B4[0, 3] = +9.0909090909090912E-002;
            B4[1, 3] = +9.0909090909090912E-002;
            B4[2, 3] = -8.1818181818181823E-001;
            B4[3, 3] = +1.0000000000000000E+000;

            Det4 = +1.6000000000000000E+002;

            I4 = MatrixBuilder.CreateMatrix(4, MatrixType.Dense);
            I4[0, 0] = +2.7500000000000002E-001;
            I4[0, 1] = -2.2500000000000001E-001;
            I4[0, 2] = +2.5000000000000001E-002;
            I4[0, 3] = +2.5000000000000001E-002;
            I4[1, 0] = +2.5000000000000001E-002;
            I4[1, 1] = +2.7500000000000002E-001;
            I4[1, 2] = -2.2500000000000001E-001;
            I4[1, 3] = +2.5000000000000001E-002;
            I4[2, 0] = +2.5000000000000001E-002;
            I4[2, 1] = +2.5000000000000001E-002;
            I4[2, 2] = +2.7500000000000002E-001;
            I4[2, 3] = -2.2500000000000001E-001;
            I4[3, 0] = -2.2500000000000001E-001;
            I4[3, 1] = +2.5000000000000001E-002;
            I4[3, 2] = +2.5000000000000001E-002;
            I4[3, 3] = +2.7500000000000002E-001;

            X4 = VectorBuilder.CreateVector(4, VectorType.Dense);
            X4[0] = +1.0000000000000000E+000;
            X4[1] = +2.0000000000000000E+000;
            X4[2] = +3.0000000000000000E+000;
            X4[3] = +4.0000000000000000E+000;

            Y4 = VectorBuilder.CreateVector(4, VectorType.Dense);
            Y4[0] = +2.0000000000000000E+001;
            Y4[1] = +2.6000000000000000E+001;
            Y4[2] = +2.8000000000000000E+001;
            Y4[3] = +2.6000000000000000E+001;


            // unit testing values - order 5

            TR5 = VectorBuilder.CreateVector(5, VectorType.Dense);
            TR5[0] = +5.0000000000000000E+000;
            TR5[1] = +4.0000000000000000E+000;
            TR5[2] = +3.0000000000000000E+000;
            TR5[3] = +2.0000000000000000E+000;
            TR5[4] = +1.0000000000000000E+000;

            LC5 = VectorBuilder.CreateVector(5, VectorType.Dense);
            LC5[0] = +5.0000000000000000E+000;
            LC5[1] = +1.0000000000000000E+000;
            LC5[2] = +2.0000000000000000E+000;
            LC5[3] = +3.0000000000000000E+000;
            LC5[4] = +4.0000000000000000E+000;

            A5 = MatrixBuilder.CreateMatrix(5, MatrixType.Dense);
            A5[0, 0] = +1.0000000000000000E+000;
            A5[1, 0] = -2.0000000000000001E-001;
            A5[1, 1] = +1.0000000000000000E+000;
            A5[2, 0] = -4.2857142857142855E-001;
            A5[2, 1] = +1.4285714285714285E-001;
            A5[2, 2] = +1.0000000000000000E+000;
            A5[3, 0] = -6.6666666666666663E-001;
            A5[3, 1] = +1.1111111111111110E-001;
            A5[3, 2] = +1.1111111111111110E-001;
            A5[3, 3] = +1.0000000000000000E+000;
            A5[4, 0] = -8.7500000000000000E-001;
            A5[4, 1] = +6.2500000000000000E-002;
            A5[4, 2] = +6.2500000000000000E-002;
            A5[4, 3] = +6.2500000000000000E-002;
            A5[4, 4] = +1.0000000000000000E+000;

            D5 = VectorBuilder.CreateVector(5, VectorType.Dense);
            D5[0] = +2.0000000000000001E-001;
            D5[1] = +2.3809523809523808E-001;
            D5[2] = +2.3333333333333334E-001;
            D5[3] = +2.2500000000000001E-001;
            D5[4] = +2.1333333333333335E-001;

            B5 = MatrixBuilder.CreateMatrix(5, MatrixType.Dense);
            B5[0, 0] = +1.0000000000000000E+000;
            B5[0, 1] = -8.0000000000000004E-001;
            B5[1, 1] = +1.0000000000000000E+000;
            B5[0, 2] = +4.7619047619047616E-002;
            B5[1, 2] = -8.0952380952380953E-001;
            B5[2, 2] = +1.0000000000000000E+000;
            B5[0, 3] = +5.5555555555555552E-002;
            B5[1, 3] = +5.5555555555555552E-002;
            B5[2, 3] = -8.3333333333333337E-001;
            B5[3, 3] = +1.0000000000000000E+000;
            B5[0, 4] = +6.2500000000000000E-002;
            B5[1, 4] = +6.2500000000000000E-002;
            B5[2, 4] = +6.2500000000000000E-002;
            B5[3, 4] = -8.7500000000000000E-001;
            B5[4, 4] = +1.0000000000000000E+000;

            Det5 = +1.8750000000000000E+003;

            I5 = MatrixBuilder.CreateMatrix(5, MatrixType.Dense);
            I5[0, 0] = +2.1333333333333335E-001;
            I5[0, 1] = -1.8666666666666668E-001;
            I5[0, 2] = +1.3333333333333334E-002;
            I5[0, 3] = +1.3333333333333334E-002;
            I5[0, 4] = +1.3333333333333334E-002;
            I5[1, 0] = +1.3333333333333334E-002;
            I5[1, 1] = +2.1333333333333335E-001;
            I5[1, 2] = -1.8666666666666668E-001;
            I5[1, 3] = +1.3333333333333334E-002;
            I5[1, 4] = +1.3333333333333334E-002;
            I5[2, 0] = +1.3333333333333334E-002;
            I5[2, 1] = +1.3333333333333334E-002;
            I5[2, 2] = +2.1333333333333335E-001;
            I5[2, 3] = -1.8666666666666668E-001;
            I5[2, 4] = +1.3333333333333334E-002;
            I5[3, 0] = +1.3333333333333334E-002;
            I5[3, 1] = +1.3333333333333334E-002;
            I5[3, 2] = +1.3333333333333334E-002;
            I5[3, 3] = +2.1333333333333335E-001;
            I5[3, 4] = -1.8666666666666668E-001;
            I5[4, 0] = -1.8666666666666668E-001;
            I5[4, 1] = +1.3333333333333334E-002;
            I5[4, 2] = +1.3333333333333334E-002;
            I5[4, 3] = +1.3333333333333334E-002;
            I5[4, 4] = +2.1333333333333335E-001;

            X5 = VectorBuilder.CreateVector(5, VectorType.Dense);
            X5[0] = +1.0000000000000000E+000;
            X5[1] = +2.0000000000000000E+000;
            X5[2] = +3.0000000000000000E+000;
            X5[3] = +4.0000000000000000E+000;
            X5[4] = +5.0000000000000000E+000;

            Y5 = VectorBuilder.CreateVector(5, VectorType.Dense);
            Y5[0] = +3.5000000000000000E+001;
            Y5[1] = +4.5000000000000000E+001;
            Y5[2] = +5.0000000000000000E+001;
            Y5[3] = +5.0000000000000000E+001;
            Y5[4] = +4.5000000000000000E+001;


            // unit testing values - order 10

            TR10 = VectorBuilder.CreateVector(10, VectorType.Dense);
            TR10[0] = +1.0000000000000000E+001;
            TR10[1] = +9.0000000000000000E+000;
            TR10[2] = +8.0000000000000000E+000;
            TR10[3] = +7.0000000000000000E+000;
            TR10[4] = +6.0000000000000000E+000;
            TR10[5] = +5.0000000000000000E+000;
            TR10[6] = +4.0000000000000000E+000;
            TR10[7] = +3.0000000000000000E+000;
            TR10[8] = +2.0000000000000000E+000;
            TR10[9] = +1.0000000000000000E+000;

            LC10 = VectorBuilder.CreateVector(10, VectorType.Dense);
            LC10[0] = +1.0000000000000000E+001;
            LC10[1] = +1.0000000000000000E+000;
            LC10[2] = +2.0000000000000000E+000;
            LC10[3] = +3.0000000000000000E+000;
            LC10[4] = +4.0000000000000000E+000;
            LC10[5] = +5.0000000000000000E+000;
            LC10[6] = +6.0000000000000000E+000;
            LC10[7] = +7.0000000000000000E+000;
            LC10[8] = +8.0000000000000000E+000;
            LC10[9] = +9.0000000000000000E+000;

            A10 = MatrixBuilder.CreateMatrix(10, MatrixType.Dense);
            A10[0, 0] = +1.0000000000000000E+000;
            A10[1, 0] = -1.0000000000000001E-001;
            A10[1, 1] = +1.0000000000000000E+000;
            A10[2, 0] = -2.0879120879120880E-001;
            A10[2, 1] = +8.7912087912087919E-002;
            A10[2, 2] = +1.0000000000000000E+000;
            A10[3, 0] = -3.2530120481927710E-001;
            A10[3, 1] = +8.4337349397590355E-002;
            A10[3, 2] = +8.4337349397590355E-002;
            A10[3, 3] = +1.0000000000000000E+000;
            A10[4, 0] = -4.4736842105263158E-001;
            A10[4, 1] = +7.8947368421052627E-002;
            A10[4, 2] = +7.8947368421052627E-002;
            A10[4, 3] = +7.8947368421052627E-002;
            A10[4, 4] = +1.0000000000000000E+000;
            A10[5, 0] = -5.7142857142857140E-001;
            A10[5, 1] = +7.1428571428571425E-002;
            A10[5, 2] = +7.1428571428571425E-002;
            A10[5, 3] = +7.1428571428571425E-002;
            A10[5, 4] = +7.1428571428571425E-002;
            A10[5, 5] = +1.0000000000000000E+000;
            A10[6, 0] = -6.9230769230769229E-001;
            A10[6, 1] = +6.1538461538461542E-002;
            A10[6, 2] = +6.1538461538461542E-002;
            A10[6, 3] = +6.1538461538461542E-002;
            A10[6, 4] = +6.1538461538461542E-002;
            A10[6, 5] = +6.1538461538461542E-002;
            A10[6, 6] = +1.0000000000000000E+000;
            A10[7, 0] = -8.0327868852459017E-001;
            A10[7, 1] = +4.9180327868852458E-002;
            A10[7, 2] = +4.9180327868852458E-002;
            A10[7, 3] = +4.9180327868852458E-002;
            A10[7, 4] = +4.9180327868852458E-002;
            A10[7, 5] = +4.9180327868852458E-002;
            A10[7, 6] = +4.9180327868852458E-002;
            A10[7, 7] = +1.0000000000000000E+000;
            A10[8, 0] = -8.9655172413793105E-001;
            A10[8, 1] = +3.4482758620689655E-002;
            A10[8, 2] = +3.4482758620689655E-002;
            A10[8, 3] = +3.4482758620689655E-002;
            A10[8, 4] = +3.4482758620689655E-002;
            A10[8, 5] = +3.4482758620689655E-002;
            A10[8, 6] = +3.4482758620689655E-002;
            A10[8, 7] = +3.4482758620689655E-002;
            A10[8, 8] = +1.0000000000000000E+000;
            A10[9, 0] = -9.6428571428571430E-001;
            A10[9, 1] = +1.7857142857142856E-002;
            A10[9, 2] = +1.7857142857142856E-002;
            A10[9, 3] = +1.7857142857142856E-002;
            A10[9, 4] = +1.7857142857142856E-002;
            A10[9, 5] = +1.7857142857142856E-002;
            A10[9, 6] = +1.7857142857142856E-002;
            A10[9, 7] = +1.7857142857142856E-002;
            A10[9, 8] = +1.7857142857142856E-002;
            A10[9, 9] = +1.0000000000000000E+000;

            D10 = VectorBuilder.CreateVector(10, VectorType.Dense);
            D10[0] = +1.0000000000000001E-001;
            D10[1] = +1.0989010989010989E-001;
            D10[2] = +1.0963855421686747E-001;
            D10[3] = +1.0921052631578948E-001;
            D10[4] = +1.0857142857142857E-001;
            D10[5] = +1.0769230769230770E-001;
            D10[6] = +1.0655737704918032E-001;
            D10[7] = +1.0517241379310345E-001;
            D10[8] = +1.0357142857142858E-001;
            D10[9] = +1.0181818181818182E-001;

            B10 = MatrixBuilder.CreateMatrix(10, MatrixType.Dense);
            B10[0, 0] = +1.0000000000000000E+000;
            B10[0, 1] = -9.0000000000000002E-001;
            B10[1, 1] = +1.0000000000000000E+000;
            B10[0, 2] = +1.0989010989010990E-002;
            B10[1, 2] = -9.0109890109890112E-001;
            B10[2, 2] = +1.0000000000000000E+000;
            B10[0, 3] = +1.2048192771084338E-002;
            B10[1, 3] = +1.2048192771084338E-002;
            B10[2, 3] = -9.0361445783132532E-001;
            B10[3, 3] = +1.0000000000000000E+000;
            B10[0, 4] = +1.3157894736842105E-002;
            B10[1, 4] = +1.3157894736842105E-002;
            B10[2, 4] = +1.3157894736842105E-002;
            B10[3, 4] = -9.0789473684210531E-001;
            B10[4, 4] = +1.0000000000000000E+000;
            B10[0, 5] = +1.4285714285714285E-002;
            B10[1, 5] = +1.4285714285714285E-002;
            B10[2, 5] = +1.4285714285714285E-002;
            B10[3, 5] = +1.4285714285714285E-002;
            B10[4, 5] = -9.1428571428571426E-001;
            B10[5, 5] = +1.0000000000000000E+000;
            B10[0, 6] = +1.5384615384615385E-002;
            B10[1, 6] = +1.5384615384615385E-002;
            B10[2, 6] = +1.5384615384615385E-002;
            B10[3, 6] = +1.5384615384615385E-002;
            B10[4, 6] = +1.5384615384615385E-002;
            B10[5, 6] = -9.2307692307692313E-001;
            B10[6, 6] = +1.0000000000000000E+000;
            B10[0, 7] = +1.6393442622950821E-002;
            B10[1, 7] = +1.6393442622950821E-002;
            B10[2, 7] = +1.6393442622950821E-002;
            B10[3, 7] = +1.6393442622950821E-002;
            B10[4, 7] = +1.6393442622950821E-002;
            B10[5, 7] = +1.6393442622950821E-002;
            B10[6, 7] = -9.3442622950819676E-001;
            B10[7, 7] = +1.0000000000000000E+000;
            B10[0, 8] = +1.7241379310344827E-002;
            B10[1, 8] = +1.7241379310344827E-002;
            B10[2, 8] = +1.7241379310344827E-002;
            B10[3, 8] = +1.7241379310344827E-002;
            B10[4, 8] = +1.7241379310344827E-002;
            B10[5, 8] = +1.7241379310344827E-002;
            B10[6, 8] = +1.7241379310344827E-002;
            B10[7, 8] = -9.4827586206896552E-001;
            B10[8, 8] = +1.0000000000000000E+000;
            B10[0, 9] = +1.7857142857142856E-002;
            B10[1, 9] = +1.7857142857142856E-002;
            B10[2, 9] = +1.7857142857142856E-002;
            B10[3, 9] = +1.7857142857142856E-002;
            B10[4, 9] = +1.7857142857142856E-002;
            B10[5, 9] = +1.7857142857142856E-002;
            B10[6, 9] = +1.7857142857142856E-002;
            B10[7, 9] = +1.7857142857142856E-002;
            B10[8, 9] = -9.6428571428571430E-001;
            B10[9, 9] = +1.0000000000000000E+000;

            Det10 = +5.5000000000000000E+009;

            I10 = MatrixBuilder.CreateMatrix(10, MatrixType.Dense);
            I10[0, 0] = +1.0181818181818182E-001;
            I10[0, 1] = -9.8181818181818176E-002;
            I10[0, 2] = +1.8181818181818182E-003;
            I10[0, 3] = +1.8181818181818182E-003;
            I10[0, 4] = +1.8181818181818182E-003;
            I10[0, 5] = +1.8181818181818182E-003;
            I10[0, 6] = +1.8181818181818182E-003;
            I10[0, 7] = +1.8181818181818182E-003;
            I10[0, 8] = +1.8181818181818182E-003;
            I10[0, 9] = +1.8181818181818182E-003;
            I10[1, 0] = +1.8181818181818182E-003;
            I10[1, 1] = +1.0181818181818182E-001;
            I10[1, 2] = -9.8181818181818176E-002;
            I10[1, 3] = +1.8181818181818182E-003;
            I10[1, 4] = +1.8181818181818182E-003;
            I10[1, 5] = +1.8181818181818182E-003;
            I10[1, 6] = +1.8181818181818182E-003;
            I10[1, 7] = +1.8181818181818182E-003;
            I10[1, 8] = +1.8181818181818182E-003;
            I10[1, 9] = +1.8181818181818182E-003;
            I10[2, 0] = +1.8181818181818182E-003;
            I10[2, 1] = +1.8181818181818182E-003;
            I10[2, 2] = +1.0181818181818182E-001;
            I10[2, 3] = -9.8181818181818176E-002;
            I10[2, 4] = +1.8181818181818182E-003;
            I10[2, 5] = +1.8181818181818182E-003;
            I10[2, 6] = +1.8181818181818182E-003;
            I10[2, 7] = +1.8181818181818182E-003;
            I10[2, 8] = +1.8181818181818182E-003;
            I10[2, 9] = +1.8181818181818182E-003;
            I10[3, 0] = +1.8181818181818182E-003;
            I10[3, 1] = +1.8181818181818182E-003;
            I10[3, 2] = +1.8181818181818182E-003;
            I10[3, 3] = +1.0181818181818182E-001;
            I10[3, 4] = -9.8181818181818176E-002;
            I10[3, 5] = +1.8181818181818182E-003;
            I10[3, 6] = +1.8181818181818182E-003;
            I10[3, 7] = +1.8181818181818182E-003;
            I10[3, 8] = +1.8181818181818182E-003;
            I10[3, 9] = +1.8181818181818182E-003;
            I10[4, 0] = +1.8181818181818182E-003;
            I10[4, 1] = +1.8181818181818182E-003;
            I10[4, 2] = +1.8181818181818182E-003;
            I10[4, 3] = +1.8181818181818182E-003;
            I10[4, 4] = +1.0181818181818182E-001;
            I10[4, 5] = -9.8181818181818176E-002;
            I10[4, 6] = +1.8181818181818182E-003;
            I10[4, 7] = +1.8181818181818182E-003;
            I10[4, 8] = +1.8181818181818182E-003;
            I10[4, 9] = +1.8181818181818182E-003;
            I10[5, 0] = +1.8181818181818182E-003;
            I10[5, 1] = +1.8181818181818182E-003;
            I10[5, 2] = +1.8181818181818182E-003;
            I10[5, 3] = +1.8181818181818182E-003;
            I10[5, 4] = +1.8181818181818182E-003;
            I10[5, 5] = +1.0181818181818182E-001;
            I10[5, 6] = -9.8181818181818176E-002;
            I10[5, 7] = +1.8181818181818182E-003;
            I10[5, 8] = +1.8181818181818182E-003;
            I10[5, 9] = +1.8181818181818182E-003;
            I10[6, 0] = +1.8181818181818182E-003;
            I10[6, 1] = +1.8181818181818182E-003;
            I10[6, 2] = +1.8181818181818182E-003;
            I10[6, 3] = +1.8181818181818182E-003;
            I10[6, 4] = +1.8181818181818182E-003;
            I10[6, 5] = +1.8181818181818182E-003;
            I10[6, 6] = +1.0181818181818182E-001;
            I10[6, 7] = -9.8181818181818176E-002;
            I10[6, 8] = +1.8181818181818182E-003;
            I10[6, 9] = +1.8181818181818182E-003;
            I10[7, 0] = +1.8181818181818182E-003;
            I10[7, 1] = +1.8181818181818182E-003;
            I10[7, 2] = +1.8181818181818182E-003;
            I10[7, 3] = +1.8181818181818182E-003;
            I10[7, 4] = +1.8181818181818182E-003;
            I10[7, 5] = +1.8181818181818182E-003;
            I10[7, 6] = +1.8181818181818182E-003;
            I10[7, 7] = +1.0181818181818182E-001;
            I10[7, 8] = -9.8181818181818176E-002;
            I10[7, 9] = +1.8181818181818182E-003;
            I10[8, 0] = +1.8181818181818182E-003;
            I10[8, 1] = +1.8181818181818182E-003;
            I10[8, 2] = +1.8181818181818182E-003;
            I10[8, 3] = +1.8181818181818182E-003;
            I10[8, 4] = +1.8181818181818182E-003;
            I10[8, 5] = +1.8181818181818182E-003;
            I10[8, 6] = +1.8181818181818182E-003;
            I10[8, 7] = +1.8181818181818182E-003;
            I10[8, 8] = +1.0181818181818182E-001;
            I10[8, 9] = -9.8181818181818176E-002;
            I10[9, 0] = -9.8181818181818176E-002;
            I10[9, 1] = +1.8181818181818182E-003;
            I10[9, 2] = +1.8181818181818182E-003;
            I10[9, 3] = +1.8181818181818182E-003;
            I10[9, 4] = +1.8181818181818182E-003;
            I10[9, 5] = +1.8181818181818182E-003;
            I10[9, 6] = +1.8181818181818182E-003;
            I10[9, 7] = +1.8181818181818182E-003;
            I10[9, 8] = +1.8181818181818182E-003;
            I10[9, 9] = +1.0181818181818182E-001;

            X10 = VectorBuilder.CreateVector(10, VectorType.Dense);
            X10[0] = +1.0000000000000000E+000;
            X10[1] = +2.0000000000000000E+000;
            X10[2] = +3.0000000000000000E+000;
            X10[3] = +4.0000000000000000E+000;
            X10[4] = +5.0000000000000000E+000;
            X10[5] = +6.0000000000000000E+000;
            X10[6] = +7.0000000000000000E+000;
            X10[7] = +8.0000000000000000E+000;
            X10[8] = +9.0000000000000000E+000;
            X10[9] = +1.0000000000000000E+001;

            Y10 = VectorBuilder.CreateVector(10, VectorType.Dense);
            Y10[0] = +2.2000000000000000E+002;
            Y10[1] = +2.6500000000000000E+002;
            Y10[2] = +3.0000000000000000E+002;
            Y10[3] = +3.2500000000000000E+002;
            Y10[4] = +3.4000000000000000E+002;
            Y10[5] = +3.4500000000000000E+002;
            Y10[6] = +3.4000000000000000E+002;
            Y10[7] = +3.2500000000000000E+002;
            Y10[8] = +3.0000000000000000E+002;
            Y10[9] = +2.6500000000000000E+002;

            // Tolerances
            Tolerance1 = 1.000E-015;
            Tolerance2 = 1.000E-015;
            Tolerance3 = 2.200E-015;
            Tolerance4 = 1.000E-014;
            Tolerance5 = 2.000E-014;
            Tolerance10 = 2.000E-013;
        }

        #endregion Test Fixture Setup

        #region Null Parameter Tests for Constructor

        // Test constructor with a null parameter
        [Test]
        [ExpectedException(typeof(System.ArgumentNullException))]
        public void NullParameterTestforConstructor1()
        {
            double[] row = null;
            double[] column = new double[1];
            Levinson dl = new Levinson(row, column);
        }

        // Test constructor with a null parameter
        [Test]
        [ExpectedException(typeof(System.ArgumentNullException))]
        public void NullParameterTestforConstructor2()
        {
            Vector row = null;
            Vector column = VectorBuilder.CreateVector(1, VectorType.Dense);
            Levinson dl = new Levinson(row, column);
        }

        // Test constructor with a null parameter
        [Test]
        [ExpectedException(typeof(System.ArgumentNullException))]
        public void NullParameterTestforConstructor3()
        {
            double[] row = new double[1];
            double[] column = null;
            Levinson dl = new Levinson(row, column);
        }

        // Test constructor with a null parameter
        [Test]
        [ExpectedException(typeof(System.ArgumentNullException))]
        public void NullParameterTestforConstructor4()
        {
            Vector row = VectorBuilder.CreateVector(1, VectorType.Dense);
            Vector column = null;
            Levinson dl = new Levinson(row, column);
        }

        #endregion  Null Parameter Tests for Constructor

        #region Zero Length Vector Tests for Constructor

        [Test]
        [ExpectedException(typeof(System.RankException))]
        public void ZeroLengthVectorTestsforConstructor1()
        {
            double[] row = new double[0];
            double[] column = new double[1];
            Levinson dl = new Levinson(row, column);
        }

        [Test]
        [ExpectedException(typeof(System.RankException))]
        public void ZeroLengthVectorTestsforConstructor2()
        {
            double[] row = new double[1];
            double[] column = new double[0];
            Levinson dl = new Levinson(row, column);
        }

        #endregion Zero Length Vector Tests for Constructor

        #region ToeplitzMatrix Member Test

        // check get matrix
        [Test]
        public void ToeplitzMatrixMember()
        {
            Levinson dl = new Levinson(LC5, TR5);
            Matrix dldm = dl.ToeplitzMatrix();

            for (int row = 0; row < TR5.Count; row++)
            {
                for (int column = 0; column < TR5.Count; column++)
                {
                    if (column < row)
                    {
                        Assert.IsTrue(dldm[row, column] == LC5[row - column]);
                    }
                    else
                    {
                        Assert.AreEqual(dldm[row, column], TR5[column - row], Constants.ACCEPTABLE_ERROR);
                    }
                }
            }
        }

        #endregion ToeplitzMatrix Member Test

        #region ToeplitzMatrix Member Result Test

        // check get matrix
        [Test]
        public void ToeplitzMatrixMemberResult()
        {
            Levinson dl = new Levinson(LC5, TR5);
            Matrix dldm = MatrixBuilder.CreateMatrix(TR5.Count, MatrixType.Dense);
            dl.ToeplitzMatrix(dldm);

            for (int row = 0; row < TR5.Count; row++)
            {
                for (int column = 0; column < TR5.Count; column++)
                {
                    if (column < row)
                    {
                        Assert.IsTrue(dldm[row, column] == LC5[row - column]);
                    }
                    else
                    {
                        Assert.AreEqual(dldm[row, column], TR5[column - row], Constants.ACCEPTABLE_ERROR);
                    }
                }
            }
        }

        #endregion ToeplitzMatrix Member Test

        #region ToeplitzMatrix Member Result Error Tests

        [Test]
        [ExpectedException(typeof(System.ArgumentNullException))]
        public void ToeplitzMatrixMemberResultNull()
        {
            Levinson dl = new Levinson(LC5, TR5);
            Matrix dldm = null;
            dl.ToeplitzMatrix(dldm);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void ToeplitzMatrixMemberResultNotConformable()
        {
            Levinson dl = new Levinson(LC5, TR5);
            Matrix dldm = MatrixBuilder.CreateMatrix(LC5.Count + 1);
            dl.ToeplitzMatrix(dldm);
        }

        #endregion ToeplitzMatrix Member Result Error Tests

        #region Diagonal Result Error Tests

        [Test]
        [ExpectedException(typeof(System.ArgumentNullException))]
        public void DiagonalResultNull()
        {
            Levinson dl = new Levinson(LC5, TR5);
            Matrix diagonal = null;
            dl.Diagonal(diagonal);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void DiagonalNotConformable()
        {
            Levinson dl = new Levinson(LC5, TR5);
            Matrix diagonal = MatrixBuilder.CreateMatrix(LC5.Count + 1);
            dl.Diagonal(diagonal);
        }

        #endregion Diagonal Result Error Tests

        #region Upper Result Error Tests

        [Test]
        [ExpectedException(typeof(System.ArgumentNullException))]
        public void UpperResultNull()
        {
            Levinson dl = new Levinson(LC5, TR5);
            Matrix upper = null;
            dl.Upper(upper);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void UpperResultNotConformable()
        {
            Levinson dl = new Levinson(LC5, TR5);
            Matrix upper = MatrixBuilder.CreateMatrix(LC5.Count + 1);
            dl.Upper(upper);
        }

        #endregion Upper Result Error Tests

        #region Lower Result Error Tests

        [Test]
        [ExpectedException(typeof(System.ArgumentNullException))]
        public void LowerResultNull()
        {
            Levinson dl = new Levinson(LC5, TR5);
            Matrix lower = null;
            dl.Lower(lower);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void LowerResultNotConformable()
        {
            Levinson dl = new Levinson(LC5, TR5);
            Matrix lower = MatrixBuilder.CreateMatrix(LC5.Count + 1);
            dl.Lower(lower);
        }

        #endregion Lower Result Error Tests

        #region Order Property Test

        // test order property
        [Test]
        public void OrderProperty()
        {
            Levinson dl = new Levinson(LC5, TR5);
            Assert.IsTrue(dl.Order == 5);
        }

        #endregion Order Property Test

        #region Decomposition Test 1

        // test the UDL factorisation for case 1
        [Test]
        public void Decomposition1()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC1, TR1);
            Matrix U = dl.Upper();
            Matrix D = dl.Diagonal();
            Matrix L = dl.Lower();

            // check the upper triangle
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    if (B1[i, j] != U[i, j])
                    {
                        e = System.Math.Abs((B1[i, j] - U[i, j]) / B1[i, j]);
                        if (e > me)
                        {
                            me = e;
                        }
                    }
                }
            }
            Assert.IsTrue(me < Tolerance1, "Maximum Error = " + me.ToString());

            // check the lower triangle
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    if (A1[i, j] != L[i, j])
                    {
                        e = System.Math.Abs((A1[i, j] - L[i, j]) / A1[i, j]);
                        if (e > me)
                        {
                            me = e;
                        }
                    }
                }
            }
            Assert.IsTrue(me < Tolerance1, "Maximum Error = " + me.ToString());

            // check the diagonal
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                e = System.Math.Abs((D1[i] - D[i, i]) / D1[i]);
                if (e > me)
                {
                    me = e;
                }
            }

            Assert.IsTrue(me < Tolerance1, "Maximum Error = " + me.ToString());
        }

        #endregion Decomposition Test 1

        #region Decomposition Test 2

        // test the UDL factorisation for case 2
        [Test]
        public void Decomposition2()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC2, TR2);
            Matrix U = dl.Upper();
            Matrix D = dl.Diagonal();
            Matrix L = dl.Lower();

            // check the upper triangle
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    if (B2[i, j] != U[i, j])
                    {
                        e = System.Math.Abs((B2[i, j] - U[i, j]) / B2[i, j]);
                        if (e > me)
                        {
                            me = e;
                        }
                    }
                }
            }
            Assert.IsTrue(me < Tolerance2, "Maximum Error = " + me.ToString());

            // check the lower triangle
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    if (A2[i, j] != L[i, j])
                    {
                        e = System.Math.Abs((A2[i, j] - L[i, j]) / A2[i, j]);
                        if (e > me)
                        {
                            me = e;
                        }
                    }
                }
            }
            Assert.IsTrue(me < Tolerance2, "Maximum Error = " + me.ToString());

            // check the diagonal
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                e = System.Math.Abs((D2[i] - D[i, i]) / D2[i]);
                if (e > me)
                {
                    me = e;
                }
            }

            Assert.IsTrue(me < Tolerance2, "Maximum Error = " + me.ToString());
        }

        #endregion Decomposition Test 2

        #region Decomposition Test 3

        // test the UDL factorisation for case 3
        [Test]
        public void Decomposition3()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC3, TR3);
            Matrix U = dl.Upper();
            Matrix D = dl.Diagonal();
            Matrix L = dl.Lower();

            // check the upper triangle
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    if (B3[i, j] != U[i, j])
                    {
                        e = System.Math.Abs((B3[i, j] - U[i, j]) / B3[i, j]);
                        if (e > me)
                        {
                            me = e;
                        }
                    }
                }
            }
            Assert.IsTrue(me < Tolerance3, "Maximum Error = " + me.ToString());

            // check the lower triangle
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    if (A3[i, j] != L[i, j])
                    {
                        e = System.Math.Abs((A3[i, j] - L[i, j]) / A3[i, j]);
                        if (e > me)
                        {
                            me = e;
                        }
                    }
                }
            }
            Assert.IsTrue(me < Tolerance3, "Maximum Error = " + me.ToString());

            // check the diagonal
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                e = System.Math.Abs((D3[i] - D[i, i]) / D3[i]);
                if (e > me)
                {
                    me = e;
                }
            }

            Assert.IsTrue(me < Tolerance3, "Maximum Error = " + me.ToString());
        }

        #endregion Decomposition Test 3

        #region Decomposition Test 4

        // test the UDL factorisation for case 4
        [Test]
        public void Decomposition4()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC4, TR4);
            Matrix U = dl.Upper();
            Matrix D = dl.Diagonal();
            Matrix L = dl.Lower();

            // check the upper triangle
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    if (B4[i, j] != U[i, j])
                    {
                        e = System.Math.Abs((B4[i, j] - U[i, j]) / B4[i, j]);
                        if (e > me)
                        {
                            me = e;
                        }
                    }
                }
            }
            Assert.IsTrue(me < Tolerance4, "Maximum Error = " + me.ToString());

            // check the lower triangle
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    if (A4[i, j] != L[i, j])
                    {
                        e = System.Math.Abs((A4[i, j] - L[i, j]) / A4[i, j]);
                        if (e > me)
                        {
                            me = e;
                        }
                    }
                }
            }
            Assert.IsTrue(me < Tolerance4, "Maximum Error = " + me.ToString());
            // check the diagonal
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                e = System.Math.Abs((D4[i] - D[i, i]) / D4[i]);
                if (e > me)
                {
                    me = e;
                }
            }

            Assert.IsTrue(me < Tolerance4, "Maximum Error = " + me.ToString());
        }

        #endregion Decomposition Test 4

        #region Decomposition Test 5

        // test the UDL factorisation for case 5
        [Test]
        public void Decomposition5()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC5, TR5);
            Matrix U = dl.Upper();
            Matrix D = dl.Diagonal();
            Matrix L = dl.Lower();

            // check the upper triangle
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    if (B5[i, j] != U[i, j])
                    {
                        e = System.Math.Abs((B5[i, j] - U[i, j]) / B5[i, j]);
                        if (e > me)
                        {
                            me = e;
                        }
                    }
                }
            }
            Assert.IsTrue(me < Tolerance5, "Maximum Error = " + me.ToString());

            // check the lower triangle
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    if (A5[i, j] != L[i, j])
                    {
                        e = System.Math.Abs((A5[i, j] - L[i, j]) / A5[i, j]);
                        if (e > me)
                        {
                            me = e;
                        }
                    }
                }
            }
            Assert.IsTrue(me < Tolerance5, "Maximum Error = " + me.ToString());
            // check the diagonal
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                e = System.Math.Abs((D5[i] - D[i, i]) / D5[i]);
                if (e > me)
                {
                    me = e;
                }
            }

            Assert.IsTrue(me < Tolerance5, "Maximum Error = " + me.ToString());
        }

        #endregion Decomposition Test 5

        #region Decomposition Test 10

        // test the UDL factorisation for case 10
        [Test]
        public void Decomposition10()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC10, TR10);
            Matrix U = dl.Upper();
            Matrix D = dl.Diagonal();
            Matrix L = dl.Lower();

            // check the upper triangle
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    if (B10[i, j] != U[i, j])
                    {
                        e = System.Math.Abs((B10[i, j] - U[i, j]) / B10[i, j]);
                        if (e > me)
                        {
                            me = e;
                        }
                    }
                }
            }
            Assert.IsTrue(me < Tolerance10, "Maximum Error = " + me.ToString());

            // check the lower triangle
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    if (A10[i, j] != L[i, j])
                    {
                        e = System.Math.Abs((A10[i, j] - L[i, j]) / A10[i, j]);
                        if (e > me)
                        {
                            me = e;
                        }
                    }
                }
            }
            Assert.IsTrue(me < Tolerance10, "Maximum Error = " + me.ToString());

            // check the diagonal
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                e = System.Math.Abs((D10[i] - D[i, i]) / D10[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance10, "Maximum Error = " + me.ToString());
        }

        #endregion Decomposition Test 10

        #region Decomposition Result Test 1

        // test the UDL factorisation for case 1
        [Test]
        public void DecompositionResult1()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC1, TR1);
            Matrix U = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Upper(U);
            Matrix D = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Diagonal(D);
            Matrix L = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Lower(L);

            // check the upper triangle
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    if (B1[i, j] != U[i, j])
                    {
                        e = System.Math.Abs((B1[i, j] - U[i, j]) / B1[i, j]);
                        if (e > me)
                        {
                            me = e;
                        }
                    }
                }
            }
            Assert.IsTrue(me < Tolerance1, "Maximum Error = " + me.ToString());

            // check the lower triangle
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    if (A1[i, j] != L[i, j])
                    {
                        e = System.Math.Abs((A1[i, j] - L[i, j]) / A1[i, j]);
                        if (e > me)
                        {
                            me = e;
                        }
                    }
                }
            }
            Assert.IsTrue(me < Tolerance1, "Maximum Error = " + me.ToString());

            // check the diagonal
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                e = System.Math.Abs((D1[i] - D[i, i]) / D1[i]);
                if (e > me)
                {
                    me = e;
                }
            }

            Assert.IsTrue(me < Tolerance1, "Maximum Error = " + me.ToString());
        }

        #endregion Decomposition Result Test 1

        #region Decomposition Result Test 2

        // test the UDL factorisation for case 2
        [Test]
        public void DecompositionResult2()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC2, TR2);
            Matrix U = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Upper(U);
            Matrix D = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Diagonal(D);
            Matrix L = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Lower(L);

            // check the upper triangle
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    if (B2[i, j] != U[i, j])
                    {
                        e = System.Math.Abs((B2[i, j] - U[i, j]) / B2[i, j]);
                        if (e > me)
                        {
                            me = e;
                        }
                    }
                }
            }
            Assert.IsTrue(me < Tolerance2, "Maximum Error = " + me.ToString());

            // check the lower triangle
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    if (A2[i, j] != L[i, j])
                    {
                        e = System.Math.Abs((A2[i, j] - L[i, j]) / A2[i, j]);
                        if (e > me)
                        {
                            me = e;
                        }
                    }
                }
            }
            Assert.IsTrue(me < Tolerance2, "Maximum Error = " + me.ToString());

            // check the diagonal
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                e = System.Math.Abs((D2[i] - D[i, i]) / D2[i]);
                if (e > me)
                {
                    me = e;
                }
            }

            Assert.IsTrue(me < Tolerance2, "Maximum Error = " + me.ToString());
        }

        #endregion Decomposition Result Test 2

        #region Decomposition Result Test 3

        // test the UDL factorisation for case 3
        [Test]
        public void DecompositionResult3()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC3, TR3);
            Matrix U = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Upper(U);
            Matrix D = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Diagonal(D);
            Matrix L = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Lower(L);

            // check the upper triangle
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    if (B3[i, j] != U[i, j])
                    {
                        e = System.Math.Abs((B3[i, j] - U[i, j]) / B3[i, j]);
                        if (e > me)
                        {
                            me = e;
                        }
                    }
                }
            }
            Assert.IsTrue(me < Tolerance3, "Maximum Error = " + me.ToString());

            // check the lower triangle
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    if (A3[i, j] != L[i, j])
                    {
                        e = System.Math.Abs((A3[i, j] - L[i, j]) / A3[i, j]);
                        if (e > me)
                        {
                            me = e;
                        }
                    }
                }
            }
            Assert.IsTrue(me < Tolerance3, "Maximum Error = " + me.ToString());

            // check the diagonal
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                e = System.Math.Abs((D3[i] - D[i, i]) / D3[i]);
                if (e > me)
                {
                    me = e;
                }
            }

            Assert.IsTrue(me < Tolerance3, "Maximum Error = " + me.ToString());
        }

        #endregion Decomposition Test 3

        #region Decomposition Result Test 4

        // test the UDL factorisation for case 4
        [Test]
        public void DecompositionResult4()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC4, TR4);
            Matrix U = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Upper(U);
            Matrix D = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Diagonal(D);
            Matrix L = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Lower(L); ;

            // check the upper triangle
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    if (B4[i, j] != U[i, j])
                    {
                        e = System.Math.Abs((B4[i, j] - U[i, j]) / B4[i, j]);
                        if (e > me)
                        {
                            me = e;
                        }
                    }
                }
            }
            Assert.IsTrue(me < Tolerance4, "Maximum Error = " + me.ToString());

            // check the lower triangle
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    if (A4[i, j] != L[i, j])
                    {
                        e = System.Math.Abs((A4[i, j] - L[i, j]) / A4[i, j]);
                        if (e > me)
                        {
                            me = e;
                        }
                    }
                }
            }
            Assert.IsTrue(me < Tolerance4, "Maximum Error = " + me.ToString());

            // check the diagonal
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                e = System.Math.Abs((D4[i] - D[i, i]) / D4[i]);
                if (e > me)
                {
                    me = e;
                }
            }

            Assert.IsTrue(me < Tolerance4, "Maximum Error = " + me.ToString());
        }

        #endregion Decomposition Result Test 4

        #region Decomposition Result Test 5

        // test the UDL factorisation for case 5
        [Test]
        public void DecompositionResult5()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC5, TR5);
            Matrix U = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Upper(U);
            Matrix D = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Diagonal(D);
            Matrix L = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Lower(L);

            // check the upper triangle
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    if (B5[i, j] != U[i, j])
                    {
                        e = System.Math.Abs((B5[i, j] - U[i, j]) / B5[i, j]);
                        if (e > me)
                        {
                            me = e;
                        }
                    }
                }
            }
            Assert.IsTrue(me < Tolerance5, "Maximum Error = " + me.ToString());

            // check the lower triangle
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    if (A5[i, j] != L[i, j])
                    {
                        e = System.Math.Abs((A5[i, j] - L[i, j]) / A5[i, j]);
                        if (e > me)
                        {
                            me = e;
                        }
                    }
                }
            }
            Assert.IsTrue(me < Tolerance5, "Maximum Error = " + me.ToString());

            // check the diagonal
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                e = System.Math.Abs((D5[i] - D[i, i]) / D5[i]);
                if (e > me)
                {
                    me = e;
                }
            }

            Assert.IsTrue(me < Tolerance5, "Maximum Error = " + me.ToString());
        }

        #endregion Decomposition Result Test 5

        #region Decomposition Result Test 10

        // test the UDL factorisation for case 10
        [Test]
        public void DecompositionResult10()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC10, TR10);
            Matrix U = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Upper(U);
            Matrix D = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Diagonal(D);
            Matrix L = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Lower(L);

            // check the upper triangle
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    if (B10[i, j] != U[i, j])
                    {
                        e = System.Math.Abs((B10[i, j] - U[i, j]) / B10[i, j]);
                        if (e > me)
                        {
                            me = e;
                        }
                    }
                }
            }
            Assert.IsTrue(me < Tolerance10, "Maximum Error = " + me.ToString());

            // check the lower triangle
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    if (A10[i, j] != L[i, j])
                    {
                        e = System.Math.Abs((A10[i, j] - L[i, j]) / A10[i, j]);
                        if (e > me)
                        {
                            me = e;
                        }
                    }
                }
            }
            Assert.IsTrue(me < Tolerance10, "Maximum Error = " + me.ToString());

            // check the diagonal
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                e = System.Math.Abs((D10[i] - D[i, i]) / D10[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance10, "Maximum Error = " + me.ToString());
        }

        #endregion Decomposition Result Test 10

        #region Singularity Test 1

        // check that singular matrix is detected
        [Test]
        public void SingularityProperty1()
        {
            Levinson dl = new Levinson(LC1, TR1);
            Assert.IsFalse(dl.IsSingular());
        }

        #endregion Singularity Test 1

        #region Singularity Test 2

        // check that singular matrix is detected
        [Test]
        public void SingularityProperty2()
        {
            Levinson dl = new Levinson(LC2, TR2);
            Assert.IsFalse(dl.IsSingular());
        }

        #endregion Singularity Test 2

        #region Singularity Test 3

        // check that singular matrix is detected
        [Test]
        public void SingularityProperty3()
        {
            Levinson dl = new Levinson(LC3, TR3);
            Assert.IsFalse(dl.IsSingular());
        }

        #endregion Singularity Test 3

        #region Singularity Test 4

        // check that singular matrix is detected
        [Test]
        public void SingularityProperty4()
        {
            Levinson dl = new Levinson(LC4, TR4);
            Assert.IsFalse(dl.IsSingular());
        }

        #endregion Singularity Test 4

        #region Singularity Test 5

        // check that singular matrix is detected
        [Test]
        public void SingularityProperty5()
        {
            Levinson dl = new Levinson(LC5, TR5);
            Assert.IsFalse(dl.IsSingular());
        }

        #endregion Singularity Test 5

        #region Singularity Test 10

        // check that singular matrix is detected
        [Test]
        public void SingularityProperty10()
        {
            Levinson dl = new Levinson(LC10, TR10);
            Assert.IsFalse(dl.IsSingular());
        }

        #endregion Singularity Test 10

        #region Singularity Test - Negative Case

        // check that singular matrix is detected
        [Test]
        public void SingularityProperty()
        {
            Vector T = VectorBuilder.CreateVector(10, VectorType.Dense);
            for (int i = 1; i < 10; i++)
            {
                T[i] = (i + 1);
            }
            T[0] = -2.0;

            Levinson dl = new Levinson(T, T);
            Assert.IsTrue(dl.IsSingular());
        }

        #endregion Singularity Test - Negative Case

        #region Determinant Method Test 1

        // Test the Determinant
        [Test]
        public void DeterminantMethod1()
        {
            // calculate determinant from diagonal
            Levinson dl = new Levinson(LC1, TR1);

            // check results match
            Double e = System.Math.Abs((dl.Determinant() - Det1) / Det1);
            Assert.IsTrue(e < Tolerance1);
        }

        #endregion Determinant Method Test 1

        #region Determinant Method Test 2

        // Test the Determinant
        [Test]
        public void DeterminantMethod2()
        {
            // calculate determinant from diagonal
            Levinson dl = new Levinson(LC2, TR2);

            // check results match
            Double e = System.Math.Abs((dl.Determinant() - Det2) / Det2);
            Assert.IsTrue(e < Tolerance2);
        }

        #endregion Determinant Method Test 2

        #region Determinant Method Test 3

        // Test the Determinant
        [Test]
        public void DeterminantMethod3()
        {
            // calculate determinant from diagonal
            Levinson dl = new Levinson(LC3, TR3);

            // check results match
            Double e = System.Math.Abs((dl.Determinant() - Det3) / Det3);
            Assert.IsTrue(e < Tolerance3);
        }

        #endregion Determinant Method Test 3

        #region Determinant Method Test 4

        // Test the Determinant
        [Test]
        public void DeterminantMethod4()
        {
            // calculate determinant from diagonal
            Levinson dl = new Levinson(LC4, TR4);

            // check results match
            Double e = System.Math.Abs((dl.Determinant() - Det4) / Det4);
            Assert.IsTrue(e < Tolerance4);
        }

        #endregion Determinant Method Test 4

        #region Determinant Method Test 5

        // Test the Determinant
        [Test]
        public void DeterminantMethod5()
        {
            // calculate determinant from diagonal
            Levinson dl = new Levinson(LC5, TR5);

            // check results match
            Double e = System.Math.Abs((dl.Determinant() - Det5) / Det5);
            Assert.IsTrue(e < Tolerance5);
        }

        #endregion Determinant Method Test 5

        #region Determinant Method Test 10

        // Test the Determinant
        [Test]
        public void DeterminantMethod10()
        {
            // calculate determinant from diagonal
            Levinson dl = new Levinson(LC10, TR10);

            // check results match
            Double e = System.Math.Abs((dl.Determinant() - Det10) / Det10);
            Assert.IsTrue(e < Tolerance10);
        }

        #endregion Determinant Method Test 10

        #region Get Inverse 1

        // calculate inverse using Inverse member
        [Test]
        public void Inverse1()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC1, TR1);

            // check inverse
            Matrix I = dl.Inverse();
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I1[i, j] - I[i, j]) / I1[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance1, "Maximum Error = " + me.ToString());
        }

        #endregion Get Inverse 1

        #region Get Inverse 2

        // calculate inverse using Inverse member
        [Test]
        public void Inverse2()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC2, TR2);

            // check inverse
            Matrix I = dl.Inverse();
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I2[i, j] - I[i, j]) / I2[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }

            Assert.IsTrue(me < Tolerance2, "Maximum Error = " + me.ToString());
        }

        #endregion Get Inverse 2

        #region Get Inverse 3

        // calculate inverse using Inverse member
        [Test]
        public void Inverse3()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC3, TR3);

            // check inverse
            Matrix I = dl.Inverse();
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I3[i, j] - I[i, j]) / I3[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance3, "Maximum Error = " + me.ToString());
        }

        #endregion Get Inverse 3

        #region Get Inverse 4

        // calculate inverse using Inverse member
        [Test]
        public void Inverse4()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC4, TR4);

            // check inverse
            Matrix I = dl.Inverse();
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I4[i, j] - I[i, j]) / I4[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance4, "Maximum Error = " + me.ToString());
        }

        #endregion Get Inverse 4

        #region Get Inverse 5

        // calculate inverse using Inverse member
        [Test]
        public void Inverse5()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC5, TR5);

            // check inverse
            Matrix I = dl.Inverse();
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I5[i, j] - I[i, j]) / I5[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance5, "Maximum Error = " + me.ToString());
        }

        #endregion Get Inverse 5

        #region Get Inverse 10

        // calculate inverse using Inverse member
        [Test]
        public void Inverse10()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC10, TR10);

            // check inverse
            Matrix I = dl.Inverse();
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I10[i, j] - I[i, j]) / I10[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance10, "Maximum Error = " + me.ToString());
        }

        #endregion Get Inverse 10

        #region Get Inverse Result 1

        // calculate inverse using Inverse member
        [Test]
        public void InverseResult1()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC1, TR1);

            // check inverse
            Matrix I = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Inverse(I);
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I1[i, j] - I[i, j]) / I1[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance1, "Maximum Error = " + me.ToString());
        }

        #endregion Get Inverse Result 1

        #region Get Inverse Result 2

        // calculate inverse using Inverse member
        [Test]
        public void InverseResult2()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC2, TR2);

            // check inverse
            Matrix I = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Inverse(I);

            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I2[i, j] - I[i, j]) / I2[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }

            Assert.IsTrue(me < Tolerance2, "Maximum Error = " + me.ToString());
        }

        #endregion Get Inverse Result 2

        #region Get Inverse Result 3

        // calculate inverse using Inverse member
        [Test]
        public void InverseResult3()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC3, TR3);

            // check inverse
            Matrix I = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Inverse(I);

            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I3[i, j] - I[i, j]) / I3[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance3, "Maximum Error = " + me.ToString());
        }

        #endregion Get Inverse Result 3

        #region Get Inverse Result 4

        // calculate inverse using Inverse member
        [Test]
        public void InverseResult4()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC4, TR4);

            // check inverse
            Matrix I = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Inverse(I);

            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I4[i, j] - I[i, j]) / I4[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance4, "Maximum Error = " + me.ToString());
        }

        #endregion Get Inverse Result 4

        #region Get Inverse Result 5

        // calculate inverse using Inverse member
        [Test]
        public void InverseResult5()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC5, TR5);

            // check inverse
            Matrix I = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Inverse(I);

            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I5[i, j] - I[i, j]) / I5[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance5, "Maximum Error = " + me.ToString());
        }

        #endregion Get Inverse Result 5

        #region Get Inverse Result 10

        // calculate inverse using Inverse member
        [Test]
        public void InverseResult10()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC10, TR10);

            // check inverse
            Matrix I = MatrixBuilder.CreateMatrix(dl.Order, MatrixType.Dense);
            dl.Inverse(I);

            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I10[i, j] - I[i, j]) / I10[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance10, "Maximum Error = " + me.ToString());
        }

        #endregion Get Inverse Result 10

        #region SolveVector 1

        // Test solving a linear system
        [Test]
        public void SolveVector1()
        {
            int i;
            double e, me;
            Levinson dl = new Levinson(LC1, TR1);
            Vector X = dl.Solve(Y1);

            // determine the maximum error
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                e = System.Math.Abs((X1[i] - X[i]) / X1[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance1, "Maximum Error = " + me.ToString());
        }

        #endregion SolveVector 1

        #region SolveVector 2

        // Test solving a linear system
        [Test]
        public void SolveVector2()
        {
            int i;
            double e, me;
            Levinson dl = new Levinson(LC2, TR2);
            Vector X = dl.Solve(Y2);

            // determine the maximum error
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                e = System.Math.Abs((X2[i] - X[i]) / X2[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance2, "Maximum Error = " + me.ToString());
        }

        #endregion SolveVector 2

        #region SolveVector 3

        // Test solving a linear system
        [Test]
        public void SolveVector3()
        {
            int i;
            double e, me;
            Levinson dl = new Levinson(LC3, TR3);
            Vector X = dl.Solve(Y3);

            // determine the maximum error
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                e = System.Math.Abs((X3[i] - X[i]) / X3[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance3, "Maximum Error = " + me.ToString());
        }

        #endregion SolveVector 3

        #region SolveVector 4

        // Test solving a linear system
        [Test]
        public void SolveVector4()
        {
            int i;
            double e, me;
            Levinson dl = new Levinson(LC4, TR4);
            Vector X = dl.Solve(Y4);

            // determine the maximum error
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                e = System.Math.Abs((X4[i] - X[i]) / X4[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance4, "Maximum Error = " + me.ToString());
        }

        #endregion SolveVector 4

        #region SolveVector 5

        // Test solving a linear system
        [Test]
        public void SolveVector5()
        {
            int i;
            double e, me;
            Levinson dl = new Levinson(LC5, TR5);
            Vector X = dl.Solve(Y5);

            // determine the maximum error
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                e = System.Math.Abs((X5[i] - X[i]) / X5[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance5, "Maximum Error = " + me.ToString());
        }

        #endregion SolveVector 5

        #region SolveVector 10

        // Test solving a linear system
        [Test]
        public void SolveVector10()
        {
            int i;
            double e, me;
            Levinson dl = new Levinson(LC10, TR10);
            Vector X = dl.Solve(Y10);

            // determine the maximum error
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                e = System.Math.Abs((X10[i] - X[i]) / X10[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance10, "Maximum Error = " + me.ToString());
        }

        #endregion SolveVector 10

        #region SolveVectorResult 1

        // Test solving a linear system
        [Test]
        public void SolveVectorResult1()
        {
            int i;
            double e, me;
            Levinson dl = new Levinson(LC1, TR1);
            Vector X = VectorBuilder.CreateVector(TR1.Count, VectorType.Dense);
            dl.Solve(Y1, X);

            // determine the maximum error
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                e = System.Math.Abs((X1[i] - X[i]) / X1[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance1, "Maximum Error = " + me.ToString());
        }

        #endregion SolveVectorResult 1

        #region SolveVectorResult 2

        // Test solving a linear system
        [Test]
        public void SolveVectorResult2()
        {
            int i;
            double e, me;
            Levinson dl = new Levinson(LC2, TR2);
            Vector X = VectorBuilder.CreateVector(TR2.Count, VectorType.Dense);
            dl.Solve(Y2, X);

            // determine the maximum error
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                e = System.Math.Abs((X2[i] - X[i]) / X2[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance2, "Maximum Error = " + me.ToString());
        }

        #endregion SolveVectorResult 2

        #region SolveVectorResult 3

        // Test solving a linear system
        [Test]
        public void SolveVectorResult3()
        {
            int i;
            double e, me;
            Levinson dl = new Levinson(LC3, TR3);
            Vector X = VectorBuilder.CreateVector(TR3.Count, VectorType.Dense);
            dl.Solve(Y3, X);

            // determine the maximum error
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                e = System.Math.Abs((X3[i] - X[i]) / X3[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance3, "Maximum Error = " + me.ToString());
        }

        #endregion SolveVectorResult 3

        #region SolveVectorResult 4

        // Test solving a linear system
        [Test]
        public void SolveVectorResult4()
        {
            int i;
            double e, me;
            Levinson dl = new Levinson(LC4, TR4);
            Vector X = VectorBuilder.CreateVector(TR4.Count, VectorType.Dense);
            dl.Solve(Y4, X);

            // determine the maximum error
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                e = System.Math.Abs((X4[i] - X[i]) / X4[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance4, "Maximum Error = " + me.ToString());
        }

        #endregion SolveVectorResult 4

        #region SolveVectorResult 5

        // Test solving a linear system
        [Test]
        public void SolveVectorResult5()
        {
            int i;
            double e, me;
            Levinson dl = new Levinson(LC5, TR5);
            Vector X = VectorBuilder.CreateVector(TR5.Count, VectorType.Dense);
            dl.Solve(Y5, X);

            // determine the maximum error
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                e = System.Math.Abs((X5[i] - X[i]) / X5[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance5, "Maximum Error = " + me.ToString());
        }

        #endregion SolveVectorResult 5

        #region SolveVectorResult 10

        // Test solving a linear system
        [Test]
        public void SolveVectorResult10()
        {
            int i;
            double e, me;
            Levinson dl = new Levinson(LC10, TR10);
            Vector X = VectorBuilder.CreateVector(TR10.Count, VectorType.Dense);
            dl.Solve(Y10, X);

            // determine the maximum error
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                e = System.Math.Abs((X10[i] - X[i]) / X10[i]);
                if (e > me)
                {
                    me = e;
                }
            }
            Assert.IsTrue(me < Tolerance10, "Maximum Error = " + me.ToString());
        }

        #endregion SolveVectorResult 10

        #region Solve Matrix 1

        // calculate inverse by solving linear equations with identity RHS
        [Test]
        public void SolveMatrix1()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC1, TR1);

            // check inverse
            Matrix I = dl.Solve(MatrixBuilder.CreateIdentityMatrix(1));
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I1[i, j] - I[i, j]) / I1[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance1, "Maximum Error = " + me.ToString());
        }

        #endregion Solve Matrix 1

        #region Solve Matrix 2

        // calculate inverse by solving linear equations with identity RHS
        [Test]
        public void SolveMatrix2()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC2, TR2);

            // check inverse
            Matrix I = dl.Solve(MatrixBuilder.CreateIdentityMatrix(2));
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I2[i, j] - I[i, j]) / I2[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance2, "Maximum Error = " + me.ToString());
        }

        #endregion Solve Matrix 2

        #region Solve Matrix 3

        // calculate inverse by solving linear equations with identity RHS
        [Test]
        public void SolveMatrix3()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC3, TR3);

            // check inverse
            Matrix I = dl.Solve(MatrixBuilder.CreateIdentityMatrix(3));
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I3[i, j] - I[i, j]) / I3[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance3, "Maximum Error = " + me.ToString());
        }

        #endregion Solve Matrix 3

        #region Solve Matrix 4

        // calculate inverse by solving linear equations with identity RHS
        [Test]
        public void SolveMatrix4()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC4, TR4);

            // check inverse
            Matrix I = dl.Solve(MatrixBuilder.CreateIdentityMatrix(4));
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I4[i, j] - I[i, j]) / I4[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance4, "Maximum Error = " + me.ToString());
        }

        #endregion Solve Matrix 4

        #region Solve Matrix 5

        // calculate inverse by solving linear equations with identity RHS
        [Test]
        public void SolveMatrix5()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC5, TR5);

            // check inverse
            Matrix I = dl.Solve(MatrixBuilder.CreateIdentityMatrix(5));
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I5[i, j] - I[i, j]) / I5[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance5, "Maximum Error = " + me.ToString());
        }

        #endregion Solve Matrix 5

        #region Solve Matrix 10

        // calculate inverse by solving linear equations with identity RHS
        [Test]
        public void SolveMatrix10()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC10, TR10);

            // check inverse
            Matrix I = dl.Solve(MatrixBuilder.CreateIdentityMatrix(10));
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I10[i, j] - I[i, j]) / I10[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance10, "Maximum Error = " + me.ToString());
        }

        #endregion Solve Matrix 10

        #region Solve Matrix Result 1

        // calculate inverse by solving linear equations with identity RHS
        [Test]
        public void SolveMatrixResult1()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC1, TR1);

            // check inverse
            Matrix I = MatrixBuilder.CreateMatrix(LC1.Count);
            dl.Solve(MatrixBuilder.CreateIdentityMatrix(1), I);
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I1[i, j] - I[i, j]) / I1[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance1, "Maximum Error = " + me.ToString());
        }

        #endregion Solve Matrix Result 1

        #region Solve Matrix Result 2

        // calculate inverse by solving linear equations with identity RHS
        [Test]
        public void SolveMatrixResult2()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC2, TR2);

            // check inverse
            Matrix I = MatrixBuilder.CreateMatrix(LC2.Count);
            dl.Solve(MatrixBuilder.CreateIdentityMatrix(2), I);
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I2[i, j] - I[i, j]) / I2[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance2, "Maximum Error = " + me.ToString());
        }

        #endregion Solve Matrix Result 2

        #region Solve Matrix Result 3

        // calculate inverse by solving linear equations with identity RHS
        [Test]
        public void SolveMatrixResult3()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC3, TR3);

            // check inverse
            Matrix I = MatrixBuilder.CreateMatrix(LC3.Count);
            dl.Solve(MatrixBuilder.CreateIdentityMatrix(3), I);
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I3[i, j] - I[i, j]) / I3[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance3, "Maximum Error = " + me.ToString());
        }

        #endregion Solve Matrix Result 3

        #region Solve Matrix Result 4

        // calculate inverse by solving linear equations with identity RHS
        [Test]
        public void SolveMatrixResult4()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC4, TR4);

            // check inverse
            Matrix I = MatrixBuilder.CreateMatrix(LC4.Count);
            dl.Solve(MatrixBuilder.CreateIdentityMatrix(4), I);
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I4[i, j] - I[i, j]) / I4[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance4, "Maximum Error = " + me.ToString());
        }

        #endregion Solve Matrix Result 4

        #region Solve Matrix Result 5

        // calculate inverse by solving linear equations with identity RHS
        [Test]
        public void SolveMatrixResult5()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC5, TR5);

            // check inverse
            Matrix I = MatrixBuilder.CreateMatrix(LC5.Count);
            dl.Solve(MatrixBuilder.CreateIdentityMatrix(5), I);
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I5[i, j] - I[i, j]) / I5[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance5, "Maximum Error = " + me.ToString());
        }

        #endregion Solve Matrix Result 5

        #region Solve Matrix Result 10

        // calculate inverse by solving linear equations with identity RHS
        [Test]
        public void SolveMatrixResult10()
        {
            int i, j;
            double e, me;
            Levinson dl = new Levinson(LC10, TR10);

            // check inverse
            Matrix I = MatrixBuilder.CreateMatrix(LC10.Count);
            dl.Solve(MatrixBuilder.CreateIdentityMatrix(10), I);
            me = 0.0;
            for (i = 0; i < dl.Order; i++)
            {
                for (j = 0; j < dl.Order; j++)
                {
                    e = System.Math.Abs((I10[i, j] - I[i, j]) / I10[i, j]);
                    if (e > me)
                    {
                        me = e;
                    }
                }
            }
            Assert.IsTrue(me < Tolerance10, "Maximum Error = " + me.ToString());
        }

        #endregion Solve Matrix Result10

        #region Solve Error Tests

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveVectorNull()
        {
            Vector vector = null;
            Levinson dl = new Levinson(LC10, TR10);
            dl.Solve(vector);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveVectorNullResult1()
        {
            Vector vector = null;
            Vector result = VectorBuilder.CreateVector(LC10.Count, VectorType.Dense);
            Levinson dl = new Levinson(LC10, TR10);
            dl.Solve(vector, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveVectorNullResult2()
        {
            Vector vector = VectorBuilder.CreateVector(LC10.Count, VectorType.Dense);
            Vector result = null;
            Levinson dl = new Levinson(LC10, TR10);
            dl.Solve(vector, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveMatrixNull()
        {
            Matrix matrix = null;
            Levinson dl = new Levinson(LC10, TR10);
            dl.Solve(matrix);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveMatrixNullResult1()
        {
            Matrix matrix = null;
            Matrix result = MatrixBuilder.CreateMatrix(LC10.Count, MatrixType.Dense);
            Levinson dl = new Levinson(LC10, TR10);
            dl.Solve(matrix, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveMatrixNullResult2()
        {
            Matrix matrix = MatrixBuilder.CreateMatrix(LC10.Count, MatrixType.Dense);
            Matrix result = null;
            Levinson dl = new Levinson(LC10, TR10);
            dl.Solve(matrix, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveVectorNotConformable()
        {
            Vector vector = VectorBuilder.CreateVector(LC10.Count + 1, VectorType.Dense);
            Levinson dl = new Levinson(LC10, TR10);
            dl.Solve(vector);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveVectorNotConformableResult1()
        {
            Vector vector = VectorBuilder.CreateVector(LC10.Count + 1, VectorType.Dense);
            Vector result = VectorBuilder.CreateVector(LC10.Count, VectorType.Dense);
            Levinson dl = new Levinson(LC10, TR10);
            dl.Solve(vector, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveVectorNotConformableResult2()
        {
            Vector vector = VectorBuilder.CreateVector(LC10.Count, VectorType.Dense);
            Vector result = VectorBuilder.CreateVector(LC10.Count + 1, VectorType.Dense);
            Levinson dl = new Levinson(LC10, TR10);
            dl.Solve(vector, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveMatrixNotConformable()
        {
            Matrix matrix = MatrixBuilder.CreateMatrix(LC10.Count + 1, MatrixType.Dense);
            Levinson dl = new Levinson(LC10, TR10);
            dl.Solve(matrix);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveMatrixNotConformableResult1()
        {
            Matrix matrix = MatrixBuilder.CreateMatrix(LC10.Count + 1, MatrixType.Dense);
            Matrix result = MatrixBuilder.CreateMatrix(LC10.Count, MatrixType.Dense);
            Levinson dl = new Levinson(LC10, TR10);
            dl.Solve(matrix, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveMatrixNotConformableResult2()
        {
            Matrix matrix = MatrixBuilder.CreateMatrix(LC10.Count, MatrixType.Dense);
            Matrix result = MatrixBuilder.CreateMatrix(LC10.Count + 1, MatrixType.Dense);
            Levinson dl = new Levinson(LC10, TR10);
            dl.Solve(matrix, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveMatrixNotConformableResult3()
        {
            Matrix matrix = MatrixBuilder.CreateMatrix(LC10.Count, MatrixType.Dense);
            Matrix result = MatrixBuilder.CreateMatrix(LC10.Count, LC10.Count + 1, MatrixType.Dense);
            Levinson dl = new Levinson(LC10, TR10);
            dl.Solve(matrix, result);
        }

        [Test]
        [ExpectedException(typeof(SingularMatrixException))]
        public void SolveVectorSingular1()
        {
            Vector T = VectorBuilder.CreateVector(10, VectorType.Dense);
            for (int i = 1; i < 10; i++)
            {
                T[i] = (i + 1);
            }
            T[0] = -2.0;

            Levinson dl = new Levinson(T, T);
            Vector vector = VectorBuilder.CreateVector(LC10.Count, VectorType.Dense);
            dl.Solve(vector);
        }

        [Test]
        [ExpectedException(typeof(SingularMatrixException))]
        public void SolveVectorSingular2()
        {
            Vector T = VectorBuilder.CreateVector(10, VectorType.Dense);
            for (int i = 1; i < 10; i++)
            {
                T[i] = (i + 1);
            }
            T[0] = -2.0;

            Levinson dl = new Levinson(T, T);
            Vector vector = VectorBuilder.CreateVector(LC10.Count, VectorType.Dense);
            Vector result = VectorBuilder.CreateVector(LC10.Count, VectorType.Dense);
            dl.Solve(vector, result);
        }

        [Test]
        [ExpectedException(typeof(SingularMatrixException))]
        public void SolveMatrixSingular1()
        {
            Vector T = VectorBuilder.CreateVector(10, VectorType.Dense);
            for (int i = 1; i < 10; i++)
            {
                T[i] = (i + 1);
            }
            T[0] = -2.0;

            Levinson dl = new Levinson(T, T);
            Matrix matrix = MatrixBuilder.CreateMatrix(LC10.Count, MatrixType.Dense);
            dl.Solve(matrix);
        }

        [Test]
        [ExpectedException(typeof(SingularMatrixException))]
        public void SolveMatrixSingular2()
        {
            Vector T = VectorBuilder.CreateVector(10, VectorType.Dense);
            for (int i = 1; i < 10; i++)
            {
                T[i] = (i + 1);
            }
            T[0] = -2.0;

            Levinson dl = new Levinson(T, T);
            Matrix matrix = MatrixBuilder.CreateMatrix(LC10.Count, MatrixType.Dense);
            Matrix result = MatrixBuilder.CreateMatrix(LC10.Count, MatrixType.Dense);
            dl.Solve(matrix, result);
        }
        #endregion
    }
}
